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We consider the modification of the Cahn-Hilliard equation when a time delay process through 
a memory function is taken into account. We then study the process of spinodal decomposition in 
fast phase transitions associated with a conserved order parameter. Finite-time memory effects are 
seen to affect the dynamics of phase transition at short times and have the effect of delaying, in a 
significant way, the process of rapid growth of the order parameter that follows a quench into the 
spinodal region. These effects are important in several systems characterized by fast processes, like 
nonequilibrium dynamics in the early universe and in relativistic heavy-ion collisions. 
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The dynamics of phase transitions depends on whether the order parameter that characterizes the different phases 
of a system is a conserved quantity or not. In different fields of physics and chemistry the dynamics of a conserved 
order parameter has been described phenomenologically by the Cahn-Hilliard (CH) equation [jj (see also Ref. (2] for 
a review). However, it is a first order equation in time and as such does not take into account memory effects (ME), 
, which may have quite important consequences on non-equilibrium dynamics for fast phase transitions. 

The CH equation is a diffusion-reaction type of equation. Diffuse processes are characterized by microscopic 
\J2 [ scattering events. The diffusion equation, in particular, describes the limiting situation where the time between the 
I/"") • scattering events is infinitesimally short and therefore does not respect causality. Since in real systems scattering events 
proceed through finite time intervals, their description through diffusion-type of equations poses serious problems in 
real physical situations, as already recognized for a very long time 3]. Another related problem with the usual diffusion 
equation is that it also leads to the breaking of the f-sum rule which gives the frequency sum of the dynamic structure 
factor 0, However, when the time scales of the ME are much smaller than any other scales, the breaking of 
causality and its effects become negligible, like in several applications with the CH equation, particularly in problems 
of metallurgy (see for example the references in [(|). 

However, memory and causality constrains cannot be ignored when the typical microscopic time scales are large 
O ,' in comparison with the other time scales characterizing the dynamics. Fast phase transitions are expected to have 
happened in the early universe and most certainly also characterize the phase transitions expected to occur in the 
i-C highly excited matter formed in relativistic heavy-ion collisions (RHIC). In the early universe such situations may have 
^ , happened when the typical microscopic time scales for relaxation, given by the inverse of the decay width associated 
k>( ' with particle dynamics, is larger than the Hubble time. This is a situation likely to be expected when describing 
GUT phase transitions or even the inflationary dynamics Q. In RHIC one expects to learn about the QCD phase 
transition 0, For instance, from the hydrodynamic analysis of the freeze-out temperature in the most central 
Au-Au collisions at 130 A GeV, the typical reaction time is around 10-20 fm/c ^(j- However, the characteristic time 
scale of the memory function in the Langevin equation which describes the dynamics near the chiral phase transition 
is predicted to be about 1 fm/c (as, for example, it is shown in Fig. 4 of Ref. Q), which is not short enough to 
be ignored. These are examples that, when analyzing e.g. the detailed dynamics related to conserved charges, may 
require the use of a related CH equation for these conserved charges that goes beyond the linear order in the time 
derivative and where ME can be accounted for. 

In order to motivate causality in the CH equation, it is useful to start our discussion by briefly reviewing the 
introduction of ME in the diffusion equation. Consider a random walk process, described by the simplest Langevin 
equation, x(t) = £(t), where £(t) represents a noise term with correlations = and (£(£)£(£')) = 2TW(t — 

if), with the intensity of the noise T. It is common to ignore ME and assume a Gaussian white noise, for which 
W(t — t') = S(t — t'). Then, the corresponding Fokker-Planck equation for the probability distribution P(x, t) is the 
usual diffusion equation. However, if one takes ME into account through a colored noise, for example of the form 
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W(t — t') = exp(— \t — t'\/j)/j, one obtains the so-called causal diffusion equation, which can be approximated by the 
partial differential equation 



d 2 d 
y—P( x ,t) + -P(x,t) 



dt 2 



at 



rv 2 p(i,t). 



(i) 



The propagation speed of the equation is defined by y/T/j and one can easily see that the propagation speed of the 
ordinary diffusion equation (7 = 0) is infinite. A causal diffusion equation like Eq. Q has been used recently [Tlj 
to discuss the evolution of conserved charges in RH1C. In the present paper we apply these ideas to dynamic phase 
transitions involving a conserved order parameter by deriving an analogous CH equation using memory functions. To 
the best of our knowledge, this is the first numerical analysis to incorporate memory constraints in a CH equation. 

Though our primary motivations are fast phase transitions in cosmology and heavy-ion collision dynamics, we here 
adopt a very general approach that can also be of relevance in other applications, like in condensed matter systems. 
We shall then consider the following general Ginzburg-Landau (GL) free energy 



F(4>) = J * 



x 



a ,_ 6,9 c 
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where is a conserved order parameter. To describe the phase transition, we set, as usual, b ~ T c — T and hence 
the GL free energy has two minima below T c . Because the order parameter is conserved it satisfies the equation of 
continuity J^</)(x,t) + V • J(x, t) = 0, where J(x, t) is a current. The ordinary CH equation follows assuming the 

irreversible current in the form J(x, t) = —TV , where T denotes here a kind of Onsager coefficient. Note that the 

irreversible current is instantaneously produced by the thermodynamic force V 5F g ^ and as such time delay effects 
are not contained in this formulation. However, the instantaneous assumption leads to results in contradiction with 
those obtained in experiments, like in the problems of heat conduction problem, spin diffusion, dielectric relaxation, 
and so on (see discussions in respect to this in Ref. 0). It is further known that sum rules are broken under the 
instantaneous assumption, as discussed in Ref. 0,O|. To overcome these difficulties, time delay or ME should be 
taken into account. 

A traditional way 0,13 to take into account ME is to generalize the current defined above as follows, 



J(x,t)=-[ dsd 3 x.'M(x - x',t - ,s)V> 
Jo 



(3) 



where M(x, t) is a memory function expressed by the correlation function of noise, as required by the fluctuation- 
dissipation theorem of second kind. Following the experience with the diffusion equation, we use a local function in 
space as M{x,t) = r e -*/^( 3 )(x)/7, where 7 is the relaxation time of the memory function. Substituting this into 
Eq. J2Jl, we obtain 



7^J(x,t) 



-j(x, t) -rv 



5F(<j>) 



(4) 



which is the analogous to the Maxwell-Cattaneo-type equation used in the heat conduction problem 0. When 
substituting this equation into the equation of continuity, we obtain the modified CH equation, 
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When 7 = 0, this equation reduces to the ordinary CH equation. Eq. JSJ) is our main result and in the following we 
examine the practical consequences of incorporating ME in the CH equation. We will also soon define precisely what 
we mean by Eq. JSJ to be causal, which will then imply a constraint condition on V and 7. 

A simple way to assess the consequences of introducing memory in the CH equation is to analyze the short-time 
dynamics of spinodal decomposition (SD), i.e., the process of phase separation following a quench into the two phase 
region of the phase diagram . Recall that under these circumstances that spinodal decomposition is characterized 
by the exponential growth of long-wavelength fluctuations in the order parameter at short times after the quench, 
leading to the formation of domains and coarsening at later times. This is in contrast to the other possible mechanism 
for phase separation, i.e. nucleation, where the process of phase transition initiates through the decay of a metastable 
state by formation of bubbles of the stable phase that grow and percolates. 

By making a linear approximation of Eq. ifKjl. valid for small amplitude initial conditions, we obtain the equation 
for the Fourier-transformed field 0(k, t), 



d 2 



d_ 

dt 



<£(k,t) fa -rk 2 (ak 2 -fo)</>(k,i). 
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Taking an initial condition with zero time derivative and 0(k, t = 0) = 0o(k), the solution of this equation can be 
written as 

t) = <Mk) e -*/^(A+e At ^ - A_e- At ^)/(2A), (7) 



where A± = 1 ± A and A = \J\ — 47rk 2 (ak 2 — b). On the other hand, the solution of the noncausal equation is 

<MM = o (k) e - rk2(Qk2 - &))t . (8) 

The sub-indexes c and nc stand for causal and noncausal. The long wavelength instability characterizing the SD 
happens for wave-numbers such that the exponentials are larger than zero. This happens for modes with momentum 
k 2 < b/a, for both the causal and noncausal solutions. For higher values of momentum, </> c (k, t) exhibits oscillation 
with relaxation, where the relaxation time is given by the momentum independent constant 27, while </> nc (k, t) shows 
only relaxation (decay), with a rate that increases infinitely with momentum. 

Similarly, we can derive an approximate solution at late times. Below the critical temperature the parameter b is 
positive and the order parameter condenses with the size y/b/c. We can then expand the order parameter around 
4>o = \/b/c, c/>(x, t) = 4>q + 0(x, t). Substituting this into the causal CH equation and ignoring the non-linear terms, 
we obtain an equation for the Fourier transformed fluctuations, </>(k, t), analogous to 10, 

a2 

7 — 0(k, t) + -0(k, t) a -rk 2 (ak 2 + 26)0(k, t), (9) 



whose solution is 



4> c (k,t) = + B^'- 1 , (10) 



where and B^ are arbitrary constants and X' ± — (—1 ± \/l~— 47rk 2 (ak 2 + 2b) ) / (27) . On the other hand, the 
solution for (f> nc (k, t), coming from the noncausal CH equation, is 

2 /l -A r< -rk 2 (2b+ak 2 )t /-, -, \ 

4> nc {k,t) = C k e ; • (11) 



From Eq. ill Oil it follows that there is a critical momentum k' c = [—6+ \Jb 2 + a/ {4rfT)]/a, such that for k = |k| < k' c , 
the relaxation time is l/(— A^_) and the long time fluctuation modes are overdamped, while for k > k' c the relaxation 
time is 27 and the relaxation is accompanied by oscillatory fluctuation modes. These features are not exhibited by 
the fluctuation mode solution coming from the noncausal CH equation, Eq. i|ll|) . which is always overdamped for 
any momentum k. Now, from the time Fourier transform of Eq. @, we obtain the expression for the frequency 
in terms of the wave-number (neglecting the complex part coming from the first order derivative in time), w(k) = 
-\/rk 2 (ak 2 + 2b)/ j. This defines the wave-number velocity for the fluctuations. Calculating it at the critical value 
k'j?, which characterizes the maximum scale in the momentum space (because the higher momentum modes have 
rapid oscillations and cancel in computing averages), the maximum wave-number velocity is found to be given by 

✓ 00 . jg 1+ /'^ n , , (12) 

where we have defined D = 2&r and the correlation length £ = \fajb. Taking the limit 7 — > in Eq. (|12fl , one has that 
w 2 (k' c ) goes to infinity, which is consistent with our previous assertion that the original CH equation violates causality. 
Taking the opposite limit of large values for 7, Eq. I|12f) gives a constraint condition relating the parameters D and 
7, 7 > D 1 for which the wave-number velocity should have limiting value one, which leads to a causal propagation. 
For values of 7 satisfying the constraint 7 > D we can always find allowed values of parameters for which Eq. (JSJ is 
causal. This explains our previous notations. 

Let us see how this also applies for instance to the spinodal modes and then use these results e.g. for the problem of 
searching a signal of a phase transition in RHIC. By considering the fastest-growing mode of the SD, which is defined 
by (using Eq. J7J|) 



d -1 + A 



dk 2 7 

we derive the time scale of the fastest mode to be 



0, (13) 



T c {k r ) = 27 (14) 
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When 7 is very small, the time scale is reduced to T nc (k r ) = 8£ 2 /D, which agrees with the time scale of the SD 
without the effect of memory • Since the modes that give the result of Eq. I|14l) are in fact the dominant spinodal 
modes and r nc < r c , this reflects itself in an overall delay of the time formation of domains, as described by the 
causal CH equation compared to the noncausal one, as the phase transition proceeds. This feature is confirmed by 
our simulations shown below. For a problem like RH1C and a possible signature of a phase transition coming from it, 
this difference in time scales can be very pronounced and lead to a striking effect that a signal, like charge fluctuations 
and domain formation, can be so much delayed that possibly could not be observed in the current experiments. For 
instance, in RHIC, the correlation length is typically £ ~ 1 fm. Using also the relation between the parameters D 
and 7 obtained for a quark plasma la, 7 = 3D, which is consistent with our constraint condition obtained from Eq. 
()12p. and considering D ~ 3.7 fm [ill Il2|. we obtain from Eq. I|f4|l that r c ~ 6.f fm, which is to be compared with 
the result r„ c ~ 2.2 fm. This represents almost a 200% difference for the time scales for the starting of the growth of 
fluctuations in Eq. (JSJ as compared to the ordinary CH equation (for 7 = 0). 

We have solved Eq. J5j numerically on a discrete spatial square lattice using a semi-implicit scheme in time, with a 
fast Fourier transform in the spatial coordinates 0] . We have checked the stability of the results by changing lattice 
spacings and time steps. In addition, for 7 > I we have also used a leap-frog algorithm and the results obtained with 
both methods agreed very well. For the noncausal 7 = equation, we used as initial condition </>(x, t — 0) a random 
distribution in space with zero average and amplitude 10~ 3 . For 7 ^ 0, we used in addition the condition that at 
t = the first-order derivative of <f> is zero. For the numerical work, Eq. © is re-parameterized to dimensionless 
variables, conveniently defined by time i = (S/T nc )t = (2b 2 T / a)t, space coordinates Xi = Xi/t;, field <fi = \J c/b(j> and 
7 = (8/r nc )7. In terms of these variables Eq. becomes function of only one parameter, 7. Eq. JHJ) was next solved 
for several values of 7. Two representatives results, for 7 = and 7 = 40.5 (for the example analyzed in the previous 
paragraph), are shown in Fig.^ 
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FIG. 1: Volume average of <j>(t) as a function of dimensionless time i for the noncausal (solid) and causal (dashed) solutions. 
The two curves that exit the top of the plot are results from the linearized theory. 

In Fig. ^ we present the time evolution of the volume average of the order parameter <f), defined as <f>(t) = 
"W Sx(^( x ' wnere N 3 is the total number of lattice points and the average is taken over different initial random 
configurations. </>(x, t) + indicates that only the positive values of the field are considered (i.e., a specific direction for 
the field has been selected). The rapid increase of (f> reflects the phenomenon of SD. The figure also shows the results 
by solving the linear equation. It shows that it performs extremely well up to and right after the spinodal growth of 
the order parameter, then justifying our previous analytical results based on the solution of the linear equation. The 
effect of 7 is seen to be more important at earlier times, consistent with the memory function used and becomes less 
important after the rapid growth of the order parameter (the spinodal explosion). It also shows the effect of a finite 
7, increasing dramatically the delay of the spinodal explosion, as predicted by our previous analytical results. The 
time for reaching equilibrium is seen to be very long, as is common with the traditional noncausal CH equation. Also 
apparent from Fig.n are the oscillations in the order parameter for a finite 7, also predicted by our previous analysis. 
This is due to the increasing importance of the second-order time derivative as compared to the first-order one as 7 
increases, i.e. as 7 increases the dissipation term becomes less important and the equation becomes more and more a 
wave-like equation. The estimated delay for the thermalization is even larger than the recent estimation [T^j for the 



5 



time delay of the relaxation of a nonconserved order parameter. 
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FIG. 2: Spherically averaged structure function as a function of dimensionless time t for a k = 0.7fc c for the noncausal (solid) 
and causal (dashed) solutions. 

We have also investigated the effect of memory for the structure function, S(k,t) = |</>(k, t)\ 2 . This quantity is 
important because it provides information on the space-time coarsening of the domains of the different phases. In 
Fig. [3 we present the results for the spherically averaged value of S for k = 0.7k c (to emphasize the fast growth 
of the long wavelength fluctuations with k < k c ). The spherical averaging for a given k was done over momenta 



k r = \/k^ + k 2 + k% such that k — 0.1A < k r < k + 0.1A, with A = 2tt/L, where L is the size of the lattice. 



Consistently with Fig. ^ this figure shows the dramatic delay for the spinodal growth. 

As conclusion, we have introduced memory into the CH equation. For a physical situation typical for the phe- 
nomenology of RHIC, we found that these ME can delay substantially the phase-separation process and consequently, 
there might not be enough time for the system to thermalize before the breakdown of the system due to expansion. 
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